{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "This tutorial gives a short overview of the AD-module included in PorePy. For an example where the AD module has been used to solve non-linear compressible flow, see the tutorial:  \"compressible_flow_with_automatic_differentiation\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sps\n",
    "\n",
    "from porepy.ad.forward_mode import Ad_array\n",
    "import porepy.ad.functions as af"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scalar AD-variables\n",
    "\n",
    "We initiate a variable $x = 2$ by giving a pair (val, jac) to the Ad_array class. val is the value at which the function will be evaluated and jac =1 since $\\frac{d x}{dx} = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Ad_array(2, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define a function $y=x^2 + 3$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = x**2 + 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To obtain the function value and the derivative we can call .val and .jac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y value is:  7\n",
      "dy/dx is:  4.0\n"
     ]
    }
   ],
   "source": [
    "print('y value is: ', y.val)\n",
    "print('dy/dx is: ', y.jac)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$y$ is also an AD variable as a function of $x$. We can use it to declare further functions, e.g., $h(x) = e^{y(x)}$. To take the exponential of an Ad_array we need to call the exponential function found in the AD module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h value is:  1096.63315843\n",
      "dh/dx is:  4386.53263371\n"
     ]
    }
   ],
   "source": [
    "h = af.exp(y)\n",
    "print('h value is: ', h.val)\n",
    "print('dh/dx is: ', h.jac)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we knew the value and jacobian of $y$ we could alternatively skip initiating $x$ and initiate $y$ directly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h value is:  1096.63315843\n",
      "dh/dx is:  4386.53263371\n"
     ]
    }
   ],
   "source": [
    "y = Ad_array(7, 4)\n",
    "h = af.exp(y)\n",
    "print('h value is: ', h.val)\n",
    "print('dh/dx is: ', h.jac)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arrays of AD-variables\n",
    "The Ad_array class also support arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Ad_array(np.array([1,2,3]), sps.diags([1,1,1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As for the scalar case, it is straight forward to define functions using normal Python programming. Let us declare the function\n",
    "$$y = Ax + x^2$$\n",
    "which has the jacobian\n",
    "$$ J(y) = A + 2 \\text{diag}(x)$$\n",
    "With this notation we mean $x^2 = [x_1^2, x_2^2, x_3^2]$, and $\\text{diag}(x)$ is a matrix with $x$ on the diagonal and zeros elsewhere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Analytic y value: \n",
      "[14 26 32]\n",
      "Analytic y jacobian:\n",
      "[[2 2 3]\n",
      " [4 4 6]\n",
      " [7 8 6]] \n",
      "\n",
      "Ad y value: \n",
      "[14 26 32]\n",
      "Ad y jacobian:\n",
      "[[ 2.  2.  3.]\n",
      " [ 4.  4.  6.]\n",
      " [ 7.  8.  6.]]\n"
     ]
    }
   ],
   "source": [
    "A = sps.csc_matrix(np.array([[0,2,3],[4,0,6],[7,8,0]]))\n",
    "y = A*x  + x**2\n",
    "\n",
    "print('Analytic y value: ')\n",
    "print(np.array([14, 26, 32]))\n",
    "print('Analytic y jacobian:')\n",
    "print(np.array([[2,2,3],[4,4,6],[7,8,6]]),'\\n')\n",
    "print('Ad y value: ')\n",
    "print(y.val)\n",
    "print('Ad y jacobian:')\n",
    "print(y.jac.A)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
